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We study a bosonic version of the Kondo lattice model with an on-site repulsion in the conduction 
band, implemented with alkali atoms in two bands of an optical lattice. Using both weak and strong- 
coupling perturbation theory, we find that at unit filling of the conduction bosons the superfiuid 
to Mott insulator transition should be accompanied by a magnetic transition from a ferromagnet 
(in the superfiuid) to a paramagnet (in the Mott insulator). Furthermore, an analytic treatment of 
Gutzwiller mean-field theory reveals that quantum spin fluctuations induced by the Kondo exchange 
cause the otherwise continuous superfiuid to Mott-insulator phase transition to be first order. We 
show that lattice separability imposes a serious constraint on proposals to exploit excited bands for 
quantum simulations, and discuss a way to overcome this constraint in the context of our model by 
using an experimentally realized non-separable lattice. A method to probe the first-order nature of 
the transition based on collapses and revivals of the matter-wave field is also discussed. 
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I. INTRODUCTION 

Ultracold atoms in optical lattices have been used to 
simulate a variety of condensed matter Hamiltonians Ij , 
with eminent successes including the simulation of both 
the Bose [2] and Fermi [3l |4] single-hand Hubbard mod- 
els. More recently, progress in controlling and stabilizing 
atoms in the excited bands of an optical lattice [SHS] has 
given rise to the exciting possibility of simulating multi- 
band condensed matter Hamiltonians, which involve a 
nontrivial interplay of spin, charge, and orbital degrees 
of freedom. These achievements have precipitated a va- 
riety of theoretical investigations into the new physics 
made possible by the orbital degrees of freedom in an 
optical lattice [SHU]- 

The Kondo lattice model (KLM), in which tightly 
bound electrons act as spinful scattering centers for elec- 
trons in a conduction band [131, a typical example of 
the type of model one would like to simulate. In the 
KLM, the orbital degree of freedom gives rise to a rich 
phase diagram that includes, e.g., magnetically ordered 
states, a heavy Fermi liquid, and unconventional super- 
conductors. In this manuscript we will revisit a version 
of the KLM first proposed in Ref. , in which the elec- 
trons are replaced by spin-^ bosons, with the spin degree 
of freedom encoded in two hyperfine states of an alkali 
atom. Our primary new finding is that, for any small but 
nonzero Kondo coupling, the typically continuous super- 
fiuid (SF) to Mott insulator (MI) phase transition be- 
comes first-order. The qualitative difference between the 
pure Hubbard and Kondo-Hubbard model, even at ar- 
bitrarily weak Kondo coupling, is reminiscent of similar 
results for the Fermi Kondo-Hubbard model obtained in 
Ref. [15|. That the inclusion of small inter-band interac- 
tions (which are often relevant in real materials) can have 
such a dramatic effect on the Bose Hubbard phase dia- 
gram underscores the importance of generalizing optical 
lattice simulations to include orbital degrees of freedom. 



The structure of the manuscript is as follows. The ex- 
perimental implementation of the Bose Kondo-Hubbard 
model is discussed in Sec. [Til This section deviates from 
the original proposal [14] in that we emphasize the essen- 
tial role of lattice non-separability in obtaining the model 
in more than one spatial dimension. We begin our inves- 
tigation of the model's phase diagram in Sec. |III[ where 
we derive effective spin Hamiltonians valid in the weak 
coupling (Sec. HI A) and strong coupling (Sec. IIIB) 



limits in order to imderstand the magnetic properties of 
the SF and MI phases. We then employ an analytic treat- 
ment of Gutzwiller mean-field theory (MFT) in Sec. IV 
to map out the ground state phase diagram. At mean- 
field level, one can observe the interplay of two compet- 
ing tendencies: Superfiuidity of the conduction bosons 
tends to spontaneously break SU(2) symmetry, whereas 
the Kondo interaction prefers a rotationally symmetric 
ground state composed of localized singlets. While the 
rigidity of the superfiuid protects it from magnetic fluc- 
tuations, we will see that these fluctuations give rise to 
a metastable MI of Kondo singlets, causing the MI to 
SF transition to become flrst-order. Such a first-order 
transition has been realized — and confirmed by Quantum 
Monte Carlo — in the spin-i boson model of Ref. [16]. In 
Sec. jV A| we consider experimental details related to dy- 
namically maintaining a two-band model. In particular 
we discuss the implementation of the model, and the rel- 
evant parameter regimes, using the non-separable lattice 
of Ref. [T7]. In Sec. VB we discuss the preparation of 



the unit-filled MI phase. We then suggest the possibility 
of experimentally observing the first-order transition by 
ramping down the lattice to enter the SF regime; due 
to the discontinuous nature of the phase transition, even 
an arbitrarily slow ramp should excite collapses and re- 
vivals of the superfiuid order parameter. In Sec. 
summarize our findings and their relevance. 
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II. THE MODEL 

Everything that follows assumes a two dimensional 
(2D) optical lattice V{r) populated by bosonic alkali 
atoms with mass m and s-wave scattering length a^. We 
assume that all atoms are in the F = 1 hyperfine mani- 
fold with mp = ±1. In ^^Rb, the similarity of scattering 
lengths for total spin F = and F — 2 collisions, to- 
gether with the ability to offset the mp = Q state via the 



quadratic Zeeman effect [71 [TH] , strongly suppresses spin 
changing collisions into the mp = state. This justifies 
considering only two internal states for times long after 
the initial preparation, and we label these two internal 
(spin) degrees of freedom by a =t, At temperatures 
sufhciently low that s-wave scattering dominates, and ne- 
glecting nearest neighbor interactions, the Hamiltonian 
describing this system is 
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In Eq. aj^ creates a spin-cr boson in a Wannier 

orbital Wa{i — Vj) of the a*'^ band, located on the lattice 
site centered at Tj. The density operator for bosons on 

site j in the a*^ band is fij^ = '^„ci\a^ja^ ^'^'^ the 
various parameters are given by 
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The notation {a, /?} ^ {7, 5} means that these two sets 
of indices cannot contain identical elements, and is intro- 
duced so that the final term in T-L contains all (and only) 
scattering processes that change the band populations on 
a given site. The hoppings Jaij are so far unspecified; we 
will only require, since we want a 2D model, that all such 
hoppings emanating from one site to its nearest neigh- 
bors are similar in magnitude. To set the overall scale 
of kinetic energy we define the effective hopping strength 

We further restrict our attention to an initial state with 
one atom per site in the lowest vibrational band {b band, 
or localized band) of the lattice, and an average of n 
atoms per site in a single excited band (a band, or con- 
duction band). The assumption that these conditions 
will be maintained dynamically, a necessary condition 
for obtaining a bosonic Kondo model, amounts to disre- 
garding all of the WajS^yS terms in Eq. (IT]). In previous 
works [51 [HI US], such an approximation nas been justi- 
fied largely by the restricted density of states in an optical 
lattice. Very roughly, the argument is that when interac- 
tions are small compared to typical band gaps, and if the 
bands themselves are narrow, then band changing col- 
lisions tend to be off resonant and strongly suppressed. 



Here we point out a notable exception to the validity 
of this reasoning, which to our knowledge has not been 
previously described in the literature. This exception af- 
fects the implementation of the Kondo-Hubbard model, 
and also of the multi-flavor models described in Ref. [H] . 




FIG. 1: (Color online) In (a) we plot an energy level diagram 
for states |00) {S band), flO) [P^ band), |01) {Py band), |20) 
{D^2 band), |02) {Dy2 band), and |11) {D^y band). The 
resonant scattering process (due to the Wa^-yS terms in H) 
that transfers atoms from the initially populated S and Dxy 
bands into the P^ and Py bands is shown in |(b)| 

For a lattice that can be written V{x,y) = Vx{x) + 
Vy{y), a Wannier orbital in the a band can be written 
Waf}{x,y) = Wa{x)wi3{y), whereas a Wannier orbital in 
the b band is woo{x,y) = wo{x)wo{y). To have a 2D 
model, it must be true that hopping in both the x and y 
directions is greater in the conduction band than in the 
localized band, and hence a, f3 > 0. Because of the sepa- 
rability of the non-interacting part of the Hamiltonian in 
the x-y basis, the state |00)(X)|a/3) [describing a single site 
with one atom in ^00(2:, y) and one atom in Wais{x, y)] is, 
even in the presence of interactions, exactly degenerate 
with the state |0/3) (g) |aO). Furthermore, these states are 
connected by the interaction term W with a matrix ele- 
ment equal to the exchange integral V between the con- 
duction and localized orbitals. The situation is depicted 
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graphically in Fig. [T] for the case where a = (3 = 1, 
and the states |00), \aO), \0f3), and \a(3) belong to the 
S, Px, Py, and Dr^y bands respectively. The Pj. and 
Py bands will be populated via collisions on a timescale 
T = 277 /V] this is unacceptable, since r will turn out to be 
the timescale for singlet formation, and one can't expect 
to see Kondo like physics if the approximations yielding 
the model break down so quickly. The way around this 
problem is to use an optical lattice potential which can- 
not be separated in cartesian coordinates. In section [VB| 
we will give an example of an existing non-separable lat- 
tice which is favorable for avoiding this problem, but for 
the time being we continue with the approximation that 
the WaiB-yS can be ignored. 

If the atoms in the b band are deep in the unit-filled MI 
regime, and we drop terms that are therefore constant, 
then "H can be reduced to a bosonic Kondo-Hubbard 
model [H] 

ijo- j 

+ "^V^Sja- Sjb- fJ^^flja- (3) 
3 3 

In Eq. ([3]) a chemical potential has been included to 
facilitate the forthcoming mean-field treatment, and the 
spin operators are defined by Sja = \ Yl,aa' '^la'^aa'^ja^ 
with T being a vector whose components are the Pauli 
matrices. In order to avoid a clutter of indices, we have 
defined V = Vat, Jij = Jija (and J = Ja), and U = Ua- 
Because V has the same sign as Qg, the Kondo inter- 
action is antiferromagnetic (AFM) for repulsive interac- 
tions, which is a manifestation of Hund's rule adapted for 
bosons: Antisymmetrization of the spin wave function for 
two identical bosons requires the antisymmetrization of 
their spatial wavefunction, thus keeping them apart and 
lowering the energy cost of repulsive interactions. 



III. EFFECTIVE SPIN MODELS 

In the presence of the Kondo term, one still expects the 
conduction bosons to undergo a MI to SF phase transi- 
tion as the ratio U/J is varied. However, the magnetic 
properties of these phases will be heavily influenced. The 
goal of this section is to study magnetic ordering in the 
SF and MI phases by deriving effective spin Hamiltonians 
that are accurate in either the weak coupling or strong 
coupling limits. The results presented are meant to re- 
inforce and complement the mean-field picture that will 
subsequently be developed. 



A. Weak coupling 

In the fermionic KLM at V = 0, the localized spins 
are decoupled from the conduction electrons, giving rise 



to a large spin degeneracy in the ground state manifold. 
The lifting of this degeneracy is the result of virtual exci- 
tations of the conduction electrons from below to above 
the Fermi-surface, which mediate long ranged couplings 
between the localized spins known as the Rudderman- 
Kittel-Kasuya-Yosida (RKKY) interaction [20H22] . In 
dimensions greater than one, the RKKY interaction is be- 
lieved to stabilize long-range magnetic order in the KLM 
ground state [5T]. The primary difference in the bosonic 
model is that the V = spin degeneracy extends to the 
conduction band, where the ferromagnetic superfiuid can 
point in any direction. Surprisingly, this additional free- 
dom actually simplifies the situation; the degeneracy will 
now be lifted at first (rather than second) order in V, 
and at this order the effective Hamiltonian can be solved 
exactly. 

Although we are not considering the presence of a 
confining potential in this manuscript, in order to pre- 
serve the greatest generality in the ensuing discussion we 
will state our results in terms of arbitrary single parti- 
cle eigenstates i^rij) and eigenvalues for a lattice plus 
trap. At V,U = 0, the ground state is formed by putting 
all conduction bosons in the single particle wavefunction 
V'o(j) and has a degeneracy of (TVq -I- 1) x 2^, where Af is 
the number of lattice sites (also the number of b atoms) , 
and Na is the number of a atoms. It is useful to de- 
fine the functions Qj^ = 'ip*{j)^^{l), in terms of which 
the spin density operator for a conduction boson in the 
single particle groundstate is 

^a-^E^°'EsW'«;.'- (4) 

jl (TO-' 

Projection of onto the degenerate groundstate man- 
ifold yields an effective weak coupling Hamiltonian that 
is first order in the Kondo coupling V 

H^^l=2VSl-Y,g%S^,. (5) 

3 

The first order energy shift due to U has been dropped 
because it does not lift the spin degeneracy. Equation ([s]) 
describes the so-called central spin model, which is famil- 
iar in the context of electron spin decoherence in quantum 
dots |23) . In the application to quantum dots, the model 
describes coupling between the spin of an electron in the 
dot and the nuclear spins of the atomic lattice in which 
the dot sits, with a coupling function determined by the 
square of the electron's wavefunction. In the present case, 
the condensate atoms are Schwinger bosons representing 
the central spin S^, which couples to the mutually non- 
interacting localized spins with a coupling function given 
by the square of the condensate wavefunction. The cen- 
tral spin model has been studied extensively in the liter- 
ature, and exact solutions have been obtained by Bethe 
ansatz [24j . However, for a translationally invariant sys- 
tem Hwc can easily be diagonalized by rewriting it in 
terms of the conserved quantities Sf,, Sa, and s, where Sa 
is the total spin quantum number of Sa — ^jai and 
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s is the total spin quantum number of the combine spin 
Sa + Sb- The ground state is formed when the local- 
ized spins align fully ferromagnetically (sf, = N/2), and 
then their total spin Sb couples as antiferromagnetically 
as possible to the condensate spin Sa {s — — 1|). 

Second-order perturbation theory in the Kondo ex- 
change generates a correction to the weak coupling 
Hamiltonian 

+ 2V'Sa-J2n,,S,b, (6) 

j 

where TZ^i = J^r^i^iji^r)^^ [see Appendix IaI for a 
derivation of Eq. m]. The first term in Eq. MTis the 
bosonic analog of the RKKY interaction for fermions, 
and has the same physical origin; a spin at site j scatters 
an a atom out of the condensate, which then re-enters 
the condensate upon scattering off the spin at site I, and 
in this way the two spins talk to each other. The second 
term, which renormalizes the spin couplings of "Hwc , re- 
flects the ability of the scattered boson to return to the 
condensate with its spin flipped. Such a term is absent 
for fermions because Pauli exclusion principle prevents a 
conduction fermion from returning to the Fermi sea with 
its spin flipped[3^. Unlike the RKKY term this correc- 
tion is local in space, a property that can be attributed 
to the cancellation of time reversed scattering processes 
involving more than one localized spin (the sign for such 
a scattering process depends on whether the spin flip oc- 
curs when the conduction atom is exiting or returning to 
the condensate). 

Although the first term has the same physical origin as 
RKKY for fermions, an important difference arrises due 
to the structure of the coupling function. In the fermionic 
case, the Fermi surface introduces a length scale kp^, 
at which the long-range coupling oscillates. For nonin- 
teracting bosons there is no such length scale, the cou- 
pling function TZ does not oscillate, and one can show 
that it is strictly positive at any finite separation for a 
translationally invariant system in the thermodynamic 
limit. Thus any groundstate of "Hwc automatically min- 

(2) 

imizes "Hwc , and the polarization of the localized spins 
persists to second order in V (note also that the renor- 
malization of the central spin coupling constants by the 

(2) 

second term in 'Hwc is toward larger positive values). 
It should be noted, however, that MTZjj diverges loga- 
rithmically with the system size in 2D, suggesting that 
the energy cannot actually be expanded in powers of V . 
Nevertheless, the existence of exclusively ferromagnetic 
(FM) terms in the first two orders of perturbation theory 
strongly suggests a FM ground state at weak coupling. 
We note that at unit filling, the true groundstate of Hwc 
is a singlet (s = 0), but nevertheless the superfiuid and 
the localized spins are aligned ferromagnetically within 
themselves (sa = Sb = "y)- That Sb = ^ means that all 



off-diagonal elements of the correlation function 

Xji = {Sjb ■ Sib) (7) 

obtain the maximum value of ^. Hence when we say 
that the n — 1 superfiuid phase is FM, we mean that fer- 
romagnetism exists independently within the superfiuid 
and within the localized spins. 

B. Strong coupling 

We will take the strong coupling limit to be defined 
by C/ ^ J, for any V, and for simplicity we will restrict 
our discussion to the case of commensurate filling in the 
conduction band. However, we will see that a similar 
limit arrises for V ^ J and any U when the a atoms are 
at unit filling. For integer n > 1 the ground state is a MI 
with n conduction bosons per site. The eigenstates on a 
single site follow from the addition of angular momenta 
Sja + Sjb = Sj, and therefore have total spin quantum 
number = (n± l)/2. Because the interaction is AFM, 
the eigenstate with lower total spin is the ground state, 
and the MI phase with density n must have total spin 
(n — l)/2 at each site. If on a single site we label the 
state by its total spin quantum number sj and the z 
projection of total spin s|, then the n filling MI at J = 
is given by 

|MU=(g)|(n-l)/2,sp, (8) 
j 

and it has energy per site En = ^n{n—l)~ ^{2 + n). For 
ji — 1 each site contains a singlet |0,0)j, spin excitations 
are gapped (Ag = 2V), and the ground state is FM. We 
also note that for unit filling the charge gap is Ac = 
U + V, and does not vanish as C/ — > 0. Hence the unit- 
filled MI exists whenever either U or y is large compared 
to J. 

States with n > 2 can be obtained from the singlet by 
repeated application of the a atom creation operators: 

|s-,s|)«(a]^)^7+^I(at^)^7-I|0,0),. (9) 

It is interesting to note that for a MI phase with n > 2 
the charge gap is given by Ac = U, and has no depen- 
dence on V. Therefore the limit V ^ J alone does not 
yield a MI in this case. Because the states \sJ,Sj) are 
degenerate (s| — — s~,...,s~), it is possible to derive 
an effective super-exchange Hamiltonian between the to- 
tal spins on neighboring sites by perturbation theory in 
the hoppings J^ . The calculation — though complicated 
by the existence of virtual excited states with two possi- 
ble total spin quantum numbers — is straightforward, and 
yields an effective Hamiltonian that only contains the to- 
tal spin operators Sj 

nsc^-9{n)J2jf-S,-S,. (10) 
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In Eq. (10), 



gin) = 



4n + 8 f V{n + 1)^ + C/(n + 2)^ 



(n + l)2 V 



Vn + U 



(11) 



is a strictly positive, density dependent coupling con- 
stant, and we have dropped an overall density dependent 
energy shift. The ground state is therefore FM with to- 
tal spin ^(n — 1), and, unlike for the n = 1 MI, there 
is no spin gap. Notice that this result for the total spin 
also holds for the weak coupling case, although the na- 
ture of the FM state is completely different in the two 
limits. While the inter-site spin correlations between the 
localized spins are maximized at weak coupling, in the 
strong-coupling limit we find the diminished correlation 
Xji = i(^irq:i)'^i which vanishes when n = 1. 



IV. MEAN FIELD THEORY 

As was shown by a numerical analysis of the Gutzwiller 
variational ansatz in Ref. [HI, in the presence of the lo- 
calized spins the a atoms continue to exhibit a MI to SF 
phase transition at commensurate filling. Here we adopt 
an alternate but equivalent description of the Gutzwiller 
variational ansatz as a site-decoupled mean-field theory 
(MFT) ^25J , and obtain an analytical description of the 
phase transition. The motivation for this mean-field 
treatment is that (1) it semi-quantitatively describes the 
MI to SF phase transition of bosons in _D > 1 and (2) the 
new features of the present model — the Kondo term — are 
entirely local, and thus are included without further ap- 
proximation. 

The starting point of Gutzwiller mean-field theory is 
the decoupling of the kinetic energy term 



{aD^ja + al(a, J " (aL)(«,.) (12) 



Assuming that the mean-field order parameter ipa- — 
(aja-) is translationally invariant, and defining = 
4:J/U, V = V/U, and A4 — n/U, we can rewrite the 
exact Hamiltonian as a sum over identical Hamilto- 
nians defined at each site (any one of which we call the 
mean- field Hamiltonian) : 

= -jj2ii>,ai+r.a^) + jj2r.^a 

a a 
+ ^naiha-l)+2VSa- Sb-Mfla- (13) 

The mean-field ground state is obtained by minimizing 
EQAiipa) = ("Hga) with respect to ipa-- The order param- 
eter describes the formation of superfluid coherence; if "00- 
is finite, then by virtue of having to point somewhere it 
describes a state that spontaneously breaks SU(2) sym- 
metry, as the superfluid should. On the other hand, 
i/'cr = is the signature of a Mott insulator phase, in 
which number fluctuations are suppressed and phase co- 
herence vanishes. Because "Hga is invariant under the 



U(l) gauge transformations aa,ipa e^^'aaTe^^'ipu, we 
are justified to choose both components of our spinor or- 
der parameter to be real. Physically, this choice amounts 
to confining any superfluid that arrises to live in the x — z 
plane, since the y component of the (vector) superfluid 
density 



E 



(14) 



vanishes when ijj* = ipa- It should be understood in what 
follows that this apparent broken symmetry is an artifact 
maintained for simplicity. 

The mean-field Hamiltonian can easily be solved nu- 
merically by truncating the single-site Hilbert space to 
contain no more than some finite number of a atoms. 
However, much insight can be gained by proceeding as 
far as possible analytically. If one makes the assump- 
tion that the SF to MI transition is continuous, then the 
mean-field phase boundary is described exactly by per- 
turbation theory in ip^-- We want to emphasize from the 
beginning that this assumption will later be shown to 
break down in certain parameter regimes. Nevertheless, 
with this assumption in mind one writes the mean-field 
energy as 

Ega = A{J, V, M) + B{J, V, M)^^ + 0{^^), (15) 

where "0^ = "0^ is the superfluid density, and the 
phase boundary is determined by the condition B — 0. 
That Eqa can be written as an even function of ip follows 
from the symmetries of "Hqa- For a Mott lobe of filling 
n, the groundstate at finite V is in general degenerate 
(see Sec. IIIB), and B must be found by calculating the 
lowest eigenvalue of the effective Hamiltonian 



En — E„ 



When n — 1 the MI ground state is unique (the sin- 
glet), and while degenerate perturbation theory is un- 
necessary it is not wrong (we simply end up finding the 
eigenvalue of a 1 x 1 matrix, which is of course its only 
entry). Some algebra reveals that the above effective 
Hamiltonian can be rewritten in the following sensible 
form 



U 



off' 



J^ci{V,M,n)p-J^C2{V,M,n)p-S, (17) 



where S — Sa + Si,, the coefhcients Ci and C2 are pre- 
sented in Appendix [Bj and the vector superfluid density 
p acts as a magnetic field. It is not hard to show that 
C2 is strictly positive |31j . and as a result minimization 
of HcS is achieved when S points along the magnetic 
field. This result makes perfect sense, since it amounts 
to self consistency in the direction of the order param- 
eter. Because the quantum number of total spin in the 
unperturbed groundstate is s = (n — l)/2, we find 



B{n) ^ J' [ ^ + c,{V,M,n) - C2iV,M,n)'^ 



(18) 
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FIG. 2: (Color online) (a) Mott lobes in the M — J plane, with J' plotted in units of Jc ~ [S + \/8)~^, which is where the 
single band mean-field transition occurs at unit filling. The red solid lines depict the boundaries of Mott lobes for V = 17/2, 
given by M™{J, 1/2). The blue dashed hues are the boundaries of Mott lobes for the single band boson Hubbard model, 
denoted M^{J')- These have been shifted up by V = 1/2 along the A^-axis in order to more easily compare their widths with 
those of the Kondo-Hubbard lobes. The red dotted line is the boundary Alf^(^,0), which does not agree with MniJ). In 
Fig. [2|b) we plot the phase boundary of the n — 1 MI in the J-V plane. The solid blue line demarks a continuos transition, 
while the red dotted line demarks a first-order transition. Fig. shows the superfiuid order parameter tp in the J-V plane, 
showing the discontinuity along the first-order portion of the phase transition. 



B{n) — determines the boundary of the n- filling Mott 
lobe, which we will denote by M^^iJ^V). The first 
three lobes are shown in Fig. l2ja) . Notice that the width 
of the n — 1 MI is increased (from 1 to 1 -I- V along the 
A4 axis) , while the widths of the higher filling lobes are 
unaffected. The reason for this goes back to the discus- 



sion in Sec. Ill B where we pointed out that the charge 
gap is given by Ac — C/(l + V) for the n = 1 MI and by 
Ac = [/ for an n > 2 MI. 

At this point we are in a position to see that something 
is wrong with the phase boundaries as they have been 
presented so far [the solid red lines in Fig. [2f^a)]. Explicit 
calculation yields the boundaries at zero Kondo coupling 



2A^^"(^, 0) = 2n-l-J±J{J -l)^ - AJn - 



4J 
n + l 
. 

whereas the mean-field phase boundaries for the single- 
band Bose Hubbard model are given by 



(^) = 2n - 1 - J ± ^{J -l)^ -AJn, (20) 

which only agree at J' = 0. But the conduction bosons 
of the Kondo-Hubbard model certainly are governed by 
the single band Bose Hubbard model at V = 0, since 
they don't talk to the localized spins. What has gone 
wrong? The problem is that Eq. ( jlQ] ) was derived 
by assuming that the mean-field ground state in the 
MI phase has s = (n — l)/2, which is true whenever 
V > 0. However, there are also the excited spin states 
with s = {n + l)/2, separated from the ground state 
manifold by A„ = V(n -I- 1), which have been ignored. 
Formally, such a procedure is correct for any finite V, so 
long as Jij} <C A„ is satisfied. As V — ?> 0, the range of va- 
lidity shrinks, until eventually at V = the perturbative 
results break down for any finite ^j. 

Once this issue is understood, it becomes clear that the 



phase transition must in fact be first-order, and the argu- 
ment is as follows. If we sit somewhere in the red shaded 
region of Fig. [2ja) , at V = the system is a superfiuid 
(since we are outside of the dashed blue line, which i s cor- 
rect at V = 0) , and so i?GA looks like it does in Fig. |3(b)[ 
Now as we turn on small but finite V, a local minimum 
of the energy must immediately develop &t ip = Q [see 
Fig. 3(c)| , since perturbation theory has a finite range 
of validity in and predicts ;B > [a Mott insulator, 
see the red dotted curve in Fig. [2ja)]. As V becomes 
larger, the range of validity for the perturbative result 
Eq. ( 19 ) increases, until eventually the ip = Q minimum 



must be the global minimum; at this point there is a first- 
order transition into a Mott insulator of Kondo singlets. 
Notice that the first-order transition is not tied to the de- 
velopment of a cubic term in the energy, as is often the 
case, but rather results from the negative contributions 
of even powers of tp beyond ijP' . By numerically solving 
the mean-field Hamiltonian we can obtain the first-order 
phase boundaries, defined by where the metastable MI 
and the superfiuid become degenerate. The boundary 
for the first Mott lobe is shown in Fig. [2j^b) , whereas the 
discontinuous jump in the order parameter across this 
boundary is shown in Fig. [2|^c). 

While the above argument is quite rigorous, it offers 
little physical insight into what is really going on; there 
is a more intuitive argument that explains the existence 
of the metastable MI phase and pinpoints the general 
features of our model that give rise to it. The superfiuid 
ground state has, in a sense to be made precise shortly, 
a certain rigidity against spin fiuctuations. Denoting 
the V = mean-field ground state \ {Sh),p), and letting 
p — p/tp'^ be a unit vector in the direction of the super- 
fluid, one might guess that turning on a small V causes 
the groundstate to be the singlet-like | — ^ p, p) — | ^p, — p) , 
gaining 3V/2 from the Kondo term. However, it can be 
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FIG. 3: (Color online) Schematic illustration of the mean-field energy, demonstrating the mechanism that drives the SF to MI 



phase transition first-order, (a) and (b) show the two forms that -Ega can take when V = 0, (a) in the MI phase and (b) 
the SF phase, with the order parameter spontaneously breaking SU(2) symmetry by choosing a particular direction 



In (c) 



we show a possible scenario for a small nonzero V. The superfluid is weakly affected (being rigid against spin fluctuationsj 
but when the order parameter is small magnetic fluctuations are allowed, reducing the energy and causing the formation of a 
metastable MI phase. 



seen that such a state actually has ^ J'tp'^ more kinetic 
energy than the mean-field ground state at V = 0. Hence 
when J'ljj'^ > V, rather than fluctuating its spin the su- 
perfluid will just anti-align with the localized spin, gain- 
ing only V/2 from the Kondo term. Notice that this 
rigidity is purely kinetic in origin, and that the mean- 
field ground state for small V agrees, both in form and 
energy, with the lowest order results of Sec. |IIIB[ On the 
other hand, when Jip'^ < V, the superfluid gains more 
energy from the Kondo term by fluctuating its spin than 
it looses in kinetic energy, and hence its energy decreases 
by 3V/2. The tendency of the mean-field energy to be 
lowered more for small ?/' than for large ip is what enables 
a local minimum at = to arise |32| . 



V. EXPERIMENTAL DETAILS 

A. The optical lattice 

As a specific example of how the approximations and 
parameters considered in this manuscript can be obtained 
in an optical lattice, we consider the 2D potential of Ref. 

m 

I{x,y) — —4Xo[cos kx + cos ky]'^ 

- Ii[2cos{2kx -2ip) + 2co8 2ky], (21) 

and a deep transverse confining lattice with poten- 
tial I±cos{kzz). Each unit cell consists of a bi- 
ased double well [Fig. |4]Ja)], and control of the 
bias allows for an adjustment of the overlap integral 
/ dxdy\wa{x,y)\'^\'Wb{x,y)\'^ and hence of V. This lattice 
offers a large parameter space to play with, and here we 
just give one example of a lattice configuration that could 

,2 

facilitate our model. Defining En = 2m\^ ' choose 
the parameters {loj^iT^p} = {0.9i?fl, 2.52i5fl;, O.Stt}, for 
which the deep well has a depth of ^ 21Ejf and the shal- 
low well has a depth of ^ lli?_R. 



X(x, X, 0) 




(a) 

Wa{x,x) Wbix,x) 



FIG. 4: (Color online), (a) Plot of l{x,x,0), with the wave 
functions Wa and wt shown schematically, (b) A contour plot 
of X{x,y,0), showing the primary hoppings. The resulting 
spectrum in the a band (c) can be fit to a tight-binding model 
with hoppings Jx, Jy, and Jxy [see (b)]. 



Using a transverse lattice with I± — AOEn and A = 
27rfcj^ = 780nm, and the scattering length of *^Rb, we 
find Ub/iJb « 2 X 10^, Ua/4:J w 10, and V/Ua w 0.1. 
Here Jf, is the nearest neighbor hopping matrix element 
for the b band. This is just inside the n = I mean-field 
Mott lobe, with V/Ua in the correct range to observe the 
first-order phase transition, and these parameters can be 
adjusted to exit the Mott insulator in the a band, or to 
change V/Ua- One consequence of this geometry is that 
the bands are not isotropic. The effect is very small for 
the b band, but for the a band, by fitting the numerically 
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calculated dispersion [Fig. [4jc)] to a tight binding model, 
we find primary hoppings Jx = Jy — 0.002£'fl and Jxy = 

Other important energy scales are the band gaps, and 
for the parameters given above the first 4 gaps (all those 
between the b band and the a band and the gap above the 
a band) are all at least lEji. This should be compared 
to the largest relevant interaction energy Ub, and for the 
parameters given above we find Ub ~ 0.28i?_R,. It should 
be noted that this comparison is an extremely conserva- 
tive metric of how the interaction energies compare to 
the band gaps. After all, Ub is the largest interaction 
energy in the model, and the band gap separating the b 
band from the band directly above it is larger than bEn- 



B. Observation of the first-order transition 

In order to probe the SF to MI phase transition, an 
ideal starting point would be a 2D {x — y plane) MI of 
spin triplet pairs in the lowest vibrational level of the 2D 
lattice with li — 0. The PM MI could then be achieved 
with high fidelity by ramping up Ti to establish an array 
of double wells in the x — y plane, and then using either 
Raman pulses [5] or the population swapping techniques 
of Refs. [SI [7] to populate the a band. Standard time-of- 
flight imaging, combined with band-mapping techniques 
[7] and spin selective imaging, could be used to resolve 
both the superfluid coherence and the magnetic ordering 
in either band. 
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FIG. 5: (Color online) Dynamics during a slow ramp of V 
from the SF to the MI at fixed conduction band filling n = 1 
(through a first-order phase transition). The ramping func- 
tion is V(t) = V/(tanh(i/to) -I- l)/2, with to « ^.7h/ J . For a 
ramp to V/ = Vc + 5 (V/ = Vc — 5), the blue solid (dotted) 
line is ip^ and the red solid (dotted) line is {Sa)- The fast 
oscillations of tjP in the non-adiabatic case occur on a time 
scale of order ^ . 



As we have discussed in Sec. V A V can be tuned 
in a double well lattice. Therefore the first-order phase 
transition should be observable by sitting just outside the 
unit-filled Mott insulator lobe and increasing V from to 
some value V/ large enough to support a Mott insulat- 
ing ground state. The discontinuous nature of the phase 
transition will cause a failure of adiabaticity for even an 
arbitrarily slow sweep of V. At the mean-field level, the 



effect of ramping V can be explored with no further ap- 
proximation by solving the time-dependent equations of 
motion that result from minimization of the Lagrangian 



C 



n 



GA/ 



(22) 



where the expectation value is taken in the single site 
Hilbert space (for practical calculations, this space must 
be truncated by cutting off the maximum number of al- 
lowed bosons) . Defining Vc to be the value of V at which 
the metastable SF solution disappears (note that this is 
not where we define the phase boundary in Sec. IV I, we 
compare a slow ramp of V from to V/ = Vc ± i5. For 
Vf < Vc, we observe a nearly adiabatic reduction of the 
SF component, whereas for V > Vc we observe collapses 
and revivals of the SF component. This is reminiscent of 
the behavior seen in Refs. [371 [5S], where a fast quench 
was studied in the single-band Bose Hubbard model, but 
here the collapses occur even for very slow lattice ramps. 
Since the collapses are pinned to a rotation of the mag- 
netization [Fig. [5], it is not surprising that they repeat 
on a time scale ~ 

Regarding the experimental feasibility of the proposed 
dynamics around the MI to SF transition, the total time 
elapsed in Fig. [5] is lOOh/J, which is comparable to the 
longest excited band decay times measured in Ref. [S] for 
the n = 1 situation. This suggests that dynamical evi- 
dence of the first-order transition, e.g. loss of adiabaticity 
or hysteresis, should be within reach of current experi- 
ments. We also note that the first-order phase transition 
could be explored by measuring local and static observ- 
ables in the trap, for instance it could manifest as a dis- 
continuity in the density profile. 



VI. SUMMARY AND CONCLUSIONS 

In real metals, it is known that when conduction elec- 
trons interact with magnetic impurities their behavior is 
drastically altered; in order to describe such systems, it 
is necessary to study many-body Hamiltonians that in- 
clude spin, charge, and orbital degrees of freedom, such 
as the KLM. The bosonic analogues of such systems are 
experimentally accessible using ultracold alkali atoms in 
non-separable optical lattices, motivating this fairly in 
depth study of the Bose Kondo-Hubbard model. 

Our primary new finding is that the SF to MI phase 
transition of the conduction bosons is qualitatively mod- 
ified by the existence of an arbitrarily small Kondo cou- 
pling to a band of localized spins. When approaching the 
MI, as the superfluid density of the conduction bosons is 
reduced, magnetic fluctuations driven by the Kondo ex- 
change induce a metastable MI of Kondo singlets, and 
as a result the associated phase transition becomes first- 
order. We expect the first-order phase transition to be 
observable in experiment, making this model a leading 
candidate for observing entirely new, strongly correlated, 
multi-band phenomena in optical lattices. 
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Appendix A: Derivation of the weak coupling 
Hamiltonian 



Here we derive the second order weak-coupling effec- 
tive Hamiltonian via perturbation theory in the Kondo 
term Hv = 2VJ2^ Sja ■ S^b 



y(2) ^ nv\n,T.){n,j:\'H 
^ So- en 



(Al) 



The state |n, S) above has one conduction boson in the 
single particle state tpn{j), the rest of the conduction 
bosons in V'o(j): ^nd a spin configuration labeled by the 
index E; only states of this form contribute at this or- 
der. Because the energy denominators have no depen- 
dence on E, the sum over S is a completeness identity in 
spin space and can be dropped. Using the basis transfor- 



mation Oj^ 
coupling as 



7t 



ripmij) we can rewrite the Kondo 



Taking advantage of common identities for pauli matri- 
ces, within the ground state manifold we can rewrite 

(qb ^ ca \ (cb qa \ ^ ^ qb ab _ ^ qa qb ^ A4^ 



which leads to 



3:1 



2v^Sa-Y.n,,s._ 



jb- 



(A5) 



Appendix B: Parameters for Heft- 

The effective Hamiltonian that yields the Mott lobe 
boundaries [Eq. [Tt] follows from tedious but straightfor- 
ward manipulation of Eq. 16 The coefficients ci and C2 
are found to be 



Cl 



1 / 71^ 2n 



1 



2n\M-V-n 1 + M-V-n 
1 



1 + M- V + Vn-n 



(Bl) 



nv^2VY,S'^- S^,Mj)Mj), (A2) and 



where 5°„„ = J2crcr' °'\na'''<ycr'a-na' ■ 1^°^ expectation values 
within the degenerate ground state manifold we have the 
equivalence 



rajl 



-1 



C2 



{n + lf 



■n? + n \M -V - n 1+M-V-n 



1 



1+ M - V + V71- n 



(B2) 
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